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I.  INTRODUCTION 

Vortices  and  vortex  wakes  have  become  a  major  theme  of  aerodynamics  research 
since  the  advent  of  the  large  aircraft  and  the  understanding  of  their  evolution  required 
an  examination  of  many  fundamental  problems  in  fluid  mechanics.  Much  of  the 
progress  made  during  the  past  two  decades  was  discussed  at  the  Symposium  on 
Aircraft  Wake  Turbulence  and  Its  Detection  [Ref.  1]  and  at  the  Aircraft  Wake  Vortices 
Conference  [Ref.  2].  Comprehensive  reviews  of  the  entire  subject  have  been  given  by 
Donaldson  and  Bilanin  [Ref.  3],  Widnall  [Ref.  4]  and  Hallock  and  Eberle  [Ref.  5]. 

These  studies,  as  well  as  numerous  others  carried  out  since  1977,  have  uncovered 
a  number  of  complex  problems  which  must  be  resolved  in  order  to  achieve  a  better 
understanding  of  the  important  features  of -trailing  vortices  in  homogeneous  and 
stratified  media.   The  principal  ones  are  as  follows  [Ref.  6]: 

(a)  Roll-up  process:  The  velocity  and  turbulence  distribution  at  anv  station 
behind  the  wins  depend  on  the  wing  section.  win2-tip  shape.  "Reynolds 
number,  wins  incidence,  and  the  distance  of  the  stationTrom  the  wing  [Ref.  7]. 
The  distributions  of  the  initial  velocity  and  turbulence,  which  influence  the 
roll-up  and  the  decay  process,  cannot  be  changed  independently.  For 
example,  a  change  in  tip  shape  changes  the  core  size,  as  well  as  the  velocity 
and  turbulence  distributions.  High  levels  of  turbulence  result  in  an  increased 
diffusion  of  vorticity,  which  in  turn  increases  the  core  size. 

(b)  Probe  sensitivity  of  the  vortices:  Flow  visualization  studies  suggest  that 
trailing  vortices  "are  extremelv  sensitive  to  disturbances  created  by  even  very 
small  probes  or  bubbles.  This  forces  one  to  use  non-intrusive  means  of 
measurements  such  as  an  LDV.  Even  then,  'vortex  wandering'  [Ref.  8],  which 
makes  the  vortices  appear  larger  than  normal  in  time-averaged  velocity 
measurements  (for  vortices  generated  by  a  wing  in  a  wind  tunnel),  or  the 
unsteadv  nature  of  the  flow  Cfor  vortices" generated  bv  a  wing  in  a  tow  basin) 
makes  the  mean  velocity  profiles  in  the  vortices  difficult  to  determine. 

(c)  Large-scale  instabilities:  The  vortices  are  seldom  observed  to  decay  awav 
owing  to  viscous  and  turbulent  dissipation,  but  are  almost  always,  "destroyed 
by  either  mutual  induction  instability  (Crow  instability  [Ref.  9J)  and/or  vortex 
breakdown.  The  Crow  instability  grows  exponentially,  and  results  either  in  a 
linking  of  the  vortex  pair  into  a  sferies  of  crude  vortex  rings  or  in  a  highly 
disorganized  intermingling  of  the  vortices.  Vortex  breakdown,  whose 
mathematical  details  have'  not  vet  been  adeauatelv  treated,  rearranges  the 
vortex  structure  and  increases  the  core  size,  turbulence,  and  energy  dissipation. 
Thus,  it  is  very  difficult  to  measure  accurately  the  trajectories  of  the  three- 
dimensional  vortices  from  their  creation  to  their  ultimate  demise. 

(d)  Reynolds  number:  Even  the  highest  Reynolds  numbers,  based  on  wing  chord, 
reached  in  wind  tunnels  or  towing  basins,  are  in  order  of  magnatude  lower 
than  what  is  possible  for  an  aircraft.  Thus,  the  scale  effects  are  not  easy  to 
assess. 

(e)  Ambient  conditions  such  as  turbulence  and  stratification  plav  major  roles  in 
the  evolution  of  vortices.  The  quantification  of  these  'effects  requires 
numerical  analysis  and  extremely  careful  experiments 


s. 
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(f)  Ground  or  free-surface  effects:  The  vortex,  pair  may  move  towards  a  risid 
boundarv  at  which,  the  no-slip  condition  must  be  satisfied  or  towards  a  free 
surface  at  which  the  zero-shear  condition  must  be  satisfied.  In  either  case,  the 
vortices  come  under  the  influence  of  their  images  and  move  accordingly. 

The  phenomenon  is  further  complicated  by  several  additional  facts.  When  the 
vortices  are  propelled  towards  a  rigid  surface,  vorticity  of  opposite  sign  is  generated  on 
the  no-slip  boundary  and  swept  towards  the  vortex  pair.  The  total  vorticity  diminishes 
very  quickly  as  vorticity  from  the  two  regions  diffuses,  the  wall  region  serving  as  a 
strong  sink  for  the  vorticity  associated  with  the  original  vortex  [Ref.  10].  The 
developement  of  a  boundary  layer  along  the  rigid  wall  may  give  rise  to  flow  separation 
for  sufficiently  high  Reynolds  numbers.  With  or  without  such  a  separation,  however, 
the  center  of  the  vortex  pair  eventually  moves  away,  or  'rebounds',  from  the  wall 
[Refs.  10-12]. 

For  the  case  of  a  zero-stress  boundary,  the  free  surface  still  acts  as  a  vorticity 
sink,  but  this  is  relatively  weak  due  to  the  absence  of  intense  oppositely-signed 
vorticity.  Thus,  in  the  absence  of  other  impeding  phenomena,  one  expects  a  mild 
interaction  between  the  vortices  and  the  free  surface  and  a  small  rebound  of  the  vortex 
center  from  the  free  surface.  However,  the  ability  of  the  free  surface  to  deform  under' 
the  influence  of  strain  fields  leads  to  a  strong  interaction  between  the  vortices  and  the 
free  surface. 

It  is  evident  from  the  foregoing  that  the  motion  and  the  life-span  of  trailing 
vortices  are  goverened  by  a  number  of  nonlinearly-dependent  complex  phenomena.  A 
number  of  experimental  and  analytical  studies  have  been  carried  out  at  the  Naval 
Postgraduate  School  by  Sarpkaya  and  his  students  [Refs.  6,13-17]  in  order  to 
investigate  the  effects  of  these  parameters  on  the  rise  and  demise  of  the  trailing  vortices 
in  homogeneous  and  density  stratified  media.  These  studies  have  clearly  identified  the 
various  demise  mechanisms  in  both  media  and  established  basic  relationships  between 
the  rise  of  vortices  and  the  governing  parameters  in  a  finite  as  well  as  effectively  infinite 
medium  [Ref.  18],  free  from  ambient  turbulence.  The  present  investigation  is  a 
continuation  of  the  foregoing  studies  and  is  confined  to  the  computer  modeling,  via  a 
finite  differencing  technique,  of  the  rise  and  demise  of  trailing  vortices  in  stratified 
media. 
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II.  EXPERIMENTAL  EQUIPMENT  AND  PROCEDURES 

A.  EQUIPMENT 

The  equipment  used  to  generate  the  trailing  vortices  has  been  extensively  used  at 
this  facility  over  the  past  three  years  [Refs.  6,14-19].  Only  the  salient  features,  most 
recent  modifications  and  the  adaptation  for  this  work  are  briefly  described  in  the 
following. 

The  system  consists  essentially  of  a  towing  basin.  The  auxiliary  components 
consist  of  the  brine  fill  system,  turbulence  management  system,  top  and  bottom 
carriages,  velocity  measuring  system,  lighting  system,  and  the  model  [Refs.  14-19].  The 
brine  fill  system  is  controlled  by  a  computer  to  provide  a  consistently  repeatable 
density  stratified  media  [Ref.  19]. 

Drains  are  provided  at  the  bottom  of  each  end  of  the  basin.  Two  parallel  rails 
are  mounted  along  the  bottom  of  the  tank.  A  carriage  rides  smoothly  on  these  rails 
and  provides  the  test  body  with  a  constant  velocity  through  the  use  of  an  endless  cable 
and  a  DC  motor.  The  velocity  of  the  model  is  measured  and  continuously  monitored 
through  the  use  of  a  magnetic  linear  displacement  transducer.  It  yields  a  signal 
proportional  to  the  velocity  of  the  model  with  an  error  of  less  than  one  percent. 

The  two  rails,  the  carriage  and  the  filling  pipes  are  located  on  or  near  the  bottom 
of  the  basin  and  under  a  turbulence  management  system.  This  system  consists  of  one 
inch  thick  polyurethane  foam  sandwiched  between  two  perforated  aluminum  plates. 

B.  MODEL 

A  rectangular  foil  (NACA  0012)  was  used  in  the  present  study.  The  foil  had  the 
following  dimensions:  B  =  4.50  in.,  bQ  =  3.53  in.,  and  c  =  2.32  in.  The  interior  of 
this  model  was  hollowed  and  used  as  dye  reservoir  to  seed  the  vortex  cores.  The  model 
was  mounted  on  its  base  by  means  of  a  thin  streamlined  aluminum  bar  with  a  cross 
section  of  a  NACA  0006  foil  and  set  at  the  desired  angle  of  attack.  As  noted  earlier. 
the  model  was  pulled  by  means  of  a  DC  motor,  pulley,  and  cable  system  at  the  desired 
speed.  Additional  details  of  the  model  construction  and  mounting  are  discussed  by 
Johnson  [Ref.  15]. 
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C       TEST  PROCEDURES 

A  model  is  placed  in  the  basin  and  the  basin  is  filled  gradually  with  a  brine 
solution  [Ref.  19]  to  the  desired  level.  In  the  cases  where  dye  was  used  to  seed  the 
vortex  pair,  the  model  was  filled  with  dye  prior  to  filling  the  basin.  After  removal  of 
trapped  air  and  after  sufficient  time  to  allow  for  the  escape  of  dissolved  air  in  the  water 
and  the  elimination  of  any  internal  currents  in  the  basin,  the  model  was  set  in  motion. 
Most  of  the  experiments  were  repeated  at  least  three  times. 

The  motion  of  the  trailing  vortices  were  recorded  on  highspeed  Polaroid  film  at 
the  test  section  (one  of  the  plexiglass  panels  near  the  middle  of  the  basin).  Each 
picture  included  two  clocks  accurate  to  0.1  seconds,  the  vertical  and  horizontal  scales 
on  the  plexiglass  window  and,  of  course,  the  side  view  of  the  trailing  vortices  as  they 
rose  from  the  model  after  formation.  The  time  interval  between  successive  pictures  is 
determined  from  the  two  clocks.  The  first  picture  always  captured  the  model  the 
instant  it  passed  through  the  test  section.  The  subsequent  pictures  (taken  at  about 
0.75  second  intervals)  captured  the  rise  and  demise  of  the  vortices.  The  vertical  rise  of 
the  vortices  is  determined  from  the  verticle  scale.  Attention  has  been  paid  to  the  fact 
that  the  vortices  are  farther  away  from  the  scale  on  the  window  and  that  the  scale 
placed  vertically  in  the  middle  of  the  test  section  does  not  exactly  correspond  to  the 
scale  marked  on  the  window  due  to  refraction  and  parallax.  The  necessary  correction 
was  made  by  photographing  a  scale  placed  in  water  in  the  middle  of  the  test  section 
together  with  a  scale  marked  on  the  window.  This  resulted  in  a  simple  conversion 
table  which  enabled  the  determination  of  the  actual  position  of  the  vortex  from  the 
scale  reading  on  the  photograph. 

The  results  are  normalized  and  plotted  in  various  forms  and  compared  with  those 
obtained  in  the  previous  runs.  Each  experiment  was  repeated  three  times  for  most  of 
the  data  presented  herein.  All  trailing  vortices  were  recorded  on  film  until  the  time 
they  have  completely  dissipated  either  due  to  aging  (diffusion  of  vorticity  due  to 
viscousity,  turbulence,  entrainment,  and  detainment),  or  due  to  Crow  instability 
(sinusoidal  instability  and  linking  of  vortices),  and,  or  due  to  vortex  breakdown  (core 
bursting).  It  was  thus-  possible  to  determine  the  life  span  of  vortex  cores  from  their 
formation  to  their  final  demise. 
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D.       COMPUTER  MODEL 

The  computer  model  (Appendix  A)  is  run  on  the  IBM  3300  main-frame  computer 
at  the  Naval  Postgraduate  School.  The  computer  code  was  written  to  run  on  either  the 
MVS  Batch  Processing  system  or  the  VM'CMS  system.  The  stratification  parameter 
specified  for  the  experimental  runs  is  used  to  determine  other  parameters  [Ref.  I7]  used 
by  the  program.  The  code  is  compiled  under  VS  Fortran  and  is  loaded  and  executed 
using  DISSPLA  version  10.0.  The  code  was  written  such  that  the  user  has  the  capacity 
to  stipulate  the  number  of  time  step  iterations  that  are  executed.  Also,  the  time  steps 
at  which  a  constant  density  contour,  complex  velocity  profile,  density  perturbation 
contour  and  the  vorticity  contour  are  plotted,  are  designated  by  the  user  prior  to 
compilation.  These  plots  are  analized  to  determine  the  normalized  rise  height  of  the 
initial  vortex  pair  until  its  demise.  These  are  then  compared  with  those  obtained 
experimentally. 

The  computer  model  is  an  intense  computational  algorithm  that  requires 
approxiamately  12.5  minutes  of  Central  Processing  Unit  (CPU)  time  per  iteration.  Due 
to  this  high  use  of  CPU  time,  the  code  has  been  designed  to  be  executed  until  stoped  at 
a  user  defined  step  iteration  and  then  save  all  necessary  data  required  to  restart  the 
model.  The  code  can  then  be  restarted  at  that  step  iteration  at  a  later  time.  The 
program  VORRS  (VORtex  ReStart)  has  been  specifically  written  to  accomplish  this 
task. 
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III.  ANALYSIS 

A.       NUMERICAL  ANALYSIS  OF  THE  VORTEX  MOTIONS 

The  generation  of  the  internal  waves  and  the  rise  and  demise  of  the  vortices  in  a 
stratified  medium  may  well  be  analyzed  through  the  use  of  the  equations  of  motion  for 
an  incompressible  fluid  for  both  laminar  and  turbulent  motions  provided  that  a  suitable 
turbulence  closure  model  is  adopted  and  the  usual  Boussinesq  approximation 
(gravitational  acceleration  is  much  larger  than  the  fluid  accelerations)  is  made.  For  the 
type  of  motions  considered  herein  the  Boussinesq  approximation  is  quite  valid  and  has 
been  used  in  the  investigation  of  all  types  of  internal  waves  in  stratified  fluids. 

For  a  two-dimensional  flow,  with  y  vertical  and  x  horizontal,  the  equations  of 
motion  are 


du  du         du  1  dp        ■  -, 

—  +u —  +  v —  =     -  • 4-  vV^u  (3.1) 

dx  dx  dy  p  dx 


and 


dv  dv  dv  1  dp  * 

—  +u  —  +  v  —  =    -g —  +  v\rv  (3.2) 

dx  dx         dy  p  dy 


and  the  equation  of  continuity 


du         dv  -    ,„  ^ 

—  +  =  0.  (3.3) 

dx        dv 


Defining  vorticitv  as  usual  bv 


^u         dv 

;;  =  _—-_  (3.4) 

dv         dx 


and  eliminating  pressure  between  Eq.  3.1  and  Eq.  3.2  ,  one  has 
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d(,         duC,       dxC,         nV        g    5p' 

—  +  —  +  —  =  vV-  C,  +  —  —  (3.o) 

ct         dx         dy  .      pn  ox 


in  which  p'  is  the  fluctuating  part  of  the  density  p  given  by 

P  -  P0  +  P(y)  +  P'(x,y,t)  (3.6) 

where  p    is  the  reference  densitv,    and  p(v)  the  initial  density  at  elevation  v  at  t  =  0. 
The  diffusion  of  densitv  is  siven  bv 


d  d  8  , 

(—  +  u—  +  v—  -  W2)p  =  0  (3.7) 

ox  ox  dy 


Combining  Eqs.  3.3.  3.6  and  3.7,  and  simplifying,  one  has 

dp'       <3up'      dvp'  dp  ,  d2p 

—  +  —-+=   -v— K+  W2p'  +  v^4  .  (3.8) 

ox.  ox  ox  ox  ox- 


It  is  convenient  to  cast  the  foregoing  equations  into  nondimensional  forms, 
scaling  each  variable  by  a  quantity  characteristic  of  its  expected  magnitude.  Two 
possible  time  scales  exist.  There  is  the  dynamic  time  scale,  which  is  the  time  a 
characteristic  length  would  be  traversed  by  a  fluid  particle  traveling  at  the  characteristic 
velocity,  and  there  is  the  buoyant  time  scale  based  on  the  natural  buoyancy  frequency 
of  the  stratified  flow,  i.e.,  the  Brunt- Vaisala  frequency  N.  as  defined  by  Daly  [Ref.  17]. 
Each  of  the  scales  gives  a  slightly  different  form  of  the  normalized  qoveming  equations. 
Dynamic  Scale: 

Introducing  U    and  L    as  the  characteristic  velocitv  and  length,  one  has 


£Lr  Ut  L'L  U 

— ,       tm  =  — ,      Re  =  -M      F    =  — ! 
K  "       K  v  (gLV 


^m  -  I-"'       lm  =  T—  "      Re  =  -7TC>      Fv  -  n*7K  ■  (3-9) 


c  c  \  =  wcy 


Normalizing  Eqs.  3.5  and  3.8 


d'C  cuC         dv   C  1      .  i     dp   ' 

di  [    oxm  5y       }  Re      'm       ¥  2  ox  {       > 
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and 


dp    '        du   p    '       <3v   p    '  ,  1      ,  d25 

Jim.  +  ( — siim  +  — m£^  =    2V     +  — ,v2p   -  +  — ]_ms  (3_ir 

di  dx  dv  m        Re'     Hm         <5v 


m  m  '  m  -  m 


in  which 


F  -  =    — £-.     Re  = 

-s-t,   n  = 

V 

N 

0 

V    . 

N2  =  A    N  2  = 

C                   T          '                 0 

g   p(y)  . 

p    dv  • 

.    _      8    ^m 

_gn2 

Equations  3.10  ans  3.11  are  valid  when  F     >  >    1,  i.e.,  when  the  bouyancy  has 
very  little  influence  on  the  nonlinear  dynamics  of  the  motion.    In  this  case,  the  density 
perturbation  acts  as  a  passive  scalar  advected  by  the  velocity  field. 
Bouyant  scale: 

Introducing  the  following  dimensionless  parameters: 


p 

u 

V 

X 

Pm  =  ' 

Po 

U        =    """■"■ 

V        —    ~— 

x     =  — 

.     u; 

c 

■       Lc 

y 

t     =  tN  , 

m              c' 

5 

*    c 

n/=- * 

(3.12) 


and  normalizing  Eqs.  3.5  and  3.S,  one  has 


dC  du   C         dv   C               F      -           dp    ' 

-m  +    p  /        m  -m  j_         m^m  \  _   v  n2r       +       '  m                                                     (3  13") 

dt  v     5x             "dv       ;         Re      ^m       ox 

m  m                  '  m                                                m 


and 


dp  du    p  (5v    p  , 

Jim    +    F  (_jn!\n    +  ^ni)   =    F  n2  (3.14) 

Ct  OX  0v 


m  m  -  m 


in  which  Fyl  Re  and  n2  are  as  defined  in  Eq.  3.1 1. 
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Equations  3.13  and  3.14  are  valid  when  Fv  <  <  1  and  buoyancy  dominates  the 
flow.  When  F  approaches  zero,  the  equations  that  result  from  Eqs-3.13  and  3.14 
describe  the  propogation  of  linear  internal  waves.  The  F  <  <  1  regime  is  of  interest 
in  the  present  investigation  because  for  submerged  bodies  of  naval  interest  F.  (the 
Froude  number)  is  about  0.01.  The  analysis  will  consider  the  full  nonlinear  equations 
given  by  Eqs.  3.13  and  3.14  rather  than  their  linearized  form  (i.e..  the  case  of  F  =  0). 
However,  the  flow  will  be  assumed  to  be  inviscid.  i.e.,  the  viscous  diffusion  will  be 
ignored  to  a  first  order  approximation.  Subsequent  analysis  will  incorporate  the  effects 
of  viscous  and  turbulent  diffusion  into  the  numerical  calculations.  The  velocities  u  and 
v  are  given  by  the  Biot-Savart  law  as 


,      ,       f  (y'  -  yK(x',y')dx'dy' 
u(x,y)  =  j -—2 (3.15) 


and 


,.(x  -  x')L(x',v')dx'dv' 
u(x,y)  =  J-- -f-r- (3.16) 


in  which  (x,y)  is  the  point  at  which  the  velocity  is  calculated  and  {x',y')  is  an  arbitrary 
point  at  which  the  vorticity  is  £(x',y'). 
The  radial  distance  r  is  given  by: 

r  =  (x-x')2  +  (y-yf.  (3.17) 

Normalizing  Eqs.  3.15  -  3.17  through  the  use  of  the  characteristic  dimensions 
given  in  Eq.  3.12 


r    (v    '  -  v    )C(x    ',v    ')dx    'dv    ' 
P  11       —    f      -m        -  m    -  ■    m    -  m         m     -  m  /-.   ir>\ 

1  vUm   ~   J  ^2 (J-18) 

m 

,•   (x™  -  ^    K(x    '.v    ')dx    'dv    ' 
F  u     =   i  m  ' m       m    ■ m  C  icn 

v    m         J  2Kr    2  ^-l     » 

m 


and 
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r   2  =  (x     -  x   ')2  +  (y     .  v   ')2  .  (3.20) 

It  is  appropriate  to  take  V  ,  the  initial  mutual  induction  velocity  of  the  vortex 
pair,  as  the  characteristic  velocity  and  the  b  ,  the  initial  separation  of  the  vorticies.  as 
the  characteristic  lensth.    Thus,  one  has: 


r                   V            V 
L    =  b    ,     u    -  V    -  — «-  ,  F    -  — °—  = *-!/,  •  (3.21) 

c        °  c         °        2-b       v        N  b        (eb  )  2 

O  CO  -     o  • 

In  dimensionless  forms,  one  has: 

r  N  b  V  t 

r     =  2_  =  2:tF    ,      n  =  ¥v  -£-2-  ,     t     =  tN    =  — 2__    .  (3.22) 

mNb  v  •  vy'm  cbF 

co  o  o    v 

in  which  N  b  /V    is  called  the  stratification  parameter  SP.    The  value  of  SP  normally 

o    o       o  r 

varies  from  zero  to  about  100.  The  larger  the  SP,  the  stronger  the  stratification.  For 
example,  for  a  model  running  at  a  speed  of  5  ft/sec  at  an  angle  of  attack,  (a)  of  10°.  a 
SP'  smaller  than  1.0  represents  weak  stratification,  a  SP  between  1  and  5  represents 
medium  stratification,  and  a  SP  larger  than  5  corresponds  to  strong  stratification. 

B.       COMPUTER  AIDED  MODELING  OF  VORTEX  MOTION 

The  computer  algorithm  for  the  numerical  solution  of  Eqs.  3.13  and  3.14  is 
provided  in  Appendix  A.  The  solution  is  done  assuming  the  case  of  inviscid  flow  and 
the  initial  vorticitv  distribution  was  assumed  to  be  Gaussian,  i.e.: 


F..  r 


2 


?m--JLTe^(-rf-j)  (-23) 

om  om 

A  computational  domain  of  5b     bv  10b     is  used  in  the  algorithm  in  order  to 

r  o         '  O  *~ 

mimick  the  physical  features  of  the  tow  basin.  The  no-slip  condition  is  satisfied  on  the 
walls  corresponding  to  the  chanel  walls.  The  zero  normal  velocity  condition  is  assumed 
at  the  free  surface  and  along  the  vertical  line  of  flow  symmetry.  All  computations  are 
done  in  the  fourth  quadrant.  The  fields  of  velocity,  vorticity.  density,  and  the 
fluctuating  components  of  the  density  are  calculated  at  each  time  step  and  plotted  at 
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regular  intervals.  The  computer  algorithm  is  quite  robust  and  can  be  applied  to  any 
stratified  and  unstratified  flow  situation.  As  the  velocity  field  is  computed  from 
analytic  expressions  (i.e.,  the  Biot-Savart  equations)  vice  an  algorithm  involving  an 
iterative  scheme,  stability  of  the  program  is  enhanced  at  a  slight  expense  in  execution 
time.  An  upwind  finite  differencing  technique  [Ref.  20]  was  used  with  an  original  grid 
dimension  of  24  by  44  nodes  with  a  spacing  of  0.25b  .  This  grid  size  was  subsequently 
increased  to  34  by  64  and  44  by  84  in  order  to  examine  the  dependency  of  the  results 
on  the  grid  size. 

The  mesh  coordinates  (x,y)  covering  the  computational  domain  {of  dimension 
m-4  by  n-4)  are: 


X  e    2,3,4 m-2} 

(3.24) 
Ye{2,3,4,...,n-2}. 


The  necessity  for  the  the  overlap  of  two  nodes  on  all  sides  of  the  computational 
domain  is  to  facilitate  satisfying  the  boundary  condiditons.  Eqs.  3.13  and  3.14  are 
solved  over  the  original  mesh  which  results  in  vorticity  and  density  information. 

The  two  complex  velocity  grids  are  derived  by  shifting  the  original  mesh  h/2  in 
the  x  direction  for  the  u  velocity  and  h/2  in  the  y  direction  for  the  v  velocity.  The 
values  of  these  new  nodes  are  determined  using  the  vorticity  derived'  for  the  original 
computational  grid  using  the  Biot-Savart  law,  Eqs.  3.15  and  3.16.  The  velocities  at  the 
shifted  grid  points  are  then  averaged  across  the  computational  grid  points  for  further 
use  in  the  finite  differencing  scheme. 


IV.  DISCUSSION  OF  RESULTS 

The  computer  code  described  in  the  previous  section  has  been  used  to  carry  out 
calculations  for  four  specific  cases.  The  first  of  these  dealt  with  the  case  of  the  rise  of 
vortices  in  a  homogeneous  unstratified  medium.  The  grid  size  as  well  as  the  core 
radius  were  varied  within  limits  to  explore  the  differences  in  the  motion  of  the  vortices. 
Figures  B.l  through  B.13  show  the  velocity  and  vorticity  contours  at  various 
normalized  times.  Clearly,  the  initial  vortex  spacing  increases  with  time  and  the 
vortices  remain  symmetrical  due  to  the  conditions  imposed  on  the  numerical  analysis. 
In  other  words,  these  vortices  cannot  be  subjected  to  sinusoidal  instability  or  core 
bursting.  Thus,  the  numerical  results  should  agree  more  closely  with  the  experimental 
cases  where  the  vortices  did  not  undergo  such  instabilities.  Figure  B.14  shows  the 
experimental  and  numerical  results.  The  experimental  data  are  separated  as  'HIGH' 
and  'LOW  because  of  the  fact  that  the  crests  and  troughs  of  the  sinusoidal  instability 
were  separately  evaluated.  The  numerical  results  for  three  normalized  core  sizes  are 
also  shown  in  Figure  B.14.  It  is  clear  that  the  core  size  has  a  pronounced  elfect  on  the 
numerical  predictions  and  that  the  results  obtained  with  rQ  =  0.08  agree  more  closely 
with  the  experimental  data.  As  noted  earlier,  the  numerical  results  should  be  compared 
with  the  average  instantaneous  height  rather  than  with  the  'HIGH'  and  'LOW  of  the 
vortices.  There  is  no  conclusive  analysis  or  data  in  the  literature  regarding  the  vortex 
core  radius.  It  has  been  shown  previously  [Ref  16]  that  the  core  radius  depends 
strongly  on  the  shape  of  the  trailing  edges  of  the  lifting  surface.  The  better  rounded 
the  edges  of  the  control  surface  the  tighter  the  winding  of  the  vortex  sheets  and  thus 
the  smaller  the  vortex  core.  It  appears  that  a  vortex  core  of  0.07  or  0.08  should  be 
used  in  future  calculations. 

Figures  B.l 5  through  B.31  show  lines  of  constant  density,  velocity,  density 
perturbations,  and  vorticity  contours  for  the  statification  parameter  of  1.0.  Fv  = 
0.02976,  r     =   0.09.  Y/b     =    -8.0  and  X.  b     =    0.5.    The  2nd  used  was  44  bv  84. 

00  o 

These  figures  show  that  at  small  times  the  vortices  rise  vertically  and  the  circulation  in 
the  flow  field  is  primarily  due  to  the  initial  vortices.  As  time  increases  (see  e.g..  Figure 
B.21b).  the  regions  of  circulation  with  countersigned  vorticity  develope  in  the  upper 
right  and  lower  left  regions  of  the  vortex.    It  is  easv  to  show  that  this  counter  vorticity 
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not  only  reduces  the  rise  velocity  of  the  original  vortex  pair  but  pushed  the  pair  against 
each  other.  Consequently,  the  countersigned  vorticity  begins  to  dominate  the  flow  (see 
e.g.,  Figure  B.2Sb)  and  the  vortex  pair  migration  stops.  With  further  increases  in  time 
the  vortex  pair  begins  to  migrate  downward. 

The  density  contours  reveal  the  same  phenomenon  in  a  different  context.  As  the 
vortices  rise,  fluid  of  greater  density  is  pushed  upwards  (see  e.g.,  Figure  B.25a)  into 
regions  of  lesser  density.  Since  such  a  migration  cannot  go  on  indefinitely,  the  vortices 
rise  to  a  maximum  height  and  then  begin  to  sink  downwards.  The  calculations  do  not 
take  into  account  the  sinusoidal  instability  and  the  vortex  breakdown.  Consequently, 
the  vortices  in  the  numerical  calculations  continue  to  exist  unimpeded  by  the  instability 
mechanisms.  In  reality,  the  vortices  begin  to  break  up  as  they  near  the  end  of  their 
maximum  migration  and  eventually  disappear. 

The  experimental  and  calculated  values  are  compared  in  Figure  B.32.  The 
correspondence  between  the  measured  and  calculated  values  is  surprisingly  good  up  to 
the  time  of  maximum  rise.  This  is  partly  because  of  the  experimental  fact  that  the  rise 
of  vortices  in  a  highly  stratified  medium  (N  b  ,'V     =    1.00,  F     =    0.02976)  is  not 

—      '  v      O    O       O  '        v  ' 

strongly  dominated  by  sinusoidal  instability  or  vortex  breakdown.  The  migration  of 
vortices  is  inhibited  primarily  due  to  the  reduction  of  vorticity  of  the  initial  vortices  and 
the  creation  of  countersigned  vorticity. 

The  numerical  calculations  were  also  applied  to  a  medium  with  a  relatively 
smaller  stratification  for  which  the  rise  of  vortices  are  closer  -to  that  of  the 
homogeneous  case  and  thus  the  effect  of  the  core  radius  is  more  pronounced.  The 
results  obtained  with  Nb  'V    =  0.50,  F    =  0.02976  and  with  r    =  0.0S  are  shown  in 

O     O         O  '        V  o 

Figures  B.33  through  B.43.  The  observations  made  earlier  regarding  the  velocity, 
vorticity,  and  density  contours  may  be  made  with  this  particular  stratification  also. 
Obviously,  the  vortices  rise  to  a  greater  height  and  the  countersigned  vorticity  is 
relatively  weaker  and  developes  at  a  later  stage.  The  calculations  with  this  particular 
stratification  have  been  repeated  by  varying  only  the  initial  core  radius.  Figure  B.44 
shows  the  representative  experimental  data  and  the  numerical  results  obtained  with 
initial  core  radii  of  0.07,  0.08.  and  0.09.  Evidently  and  as  anticipated  on  the  basis  of 
the  previous  calculations  the  initial  core  radius  of  0.07  yields  results  which  are  in 
excellent  agreement  with  those  obtained  experimental}'.  Figure  B.44  also  shows  that 
the  rise  and  demise  of  vortices  must  be  a  strong  function  of  their  initial  core  radii. 
Experiments  in  both  the  laboratory  and  the  field  have  shown  that  to  be  true  and 
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pointed  out  the  fact  that  the   life  of  a  vortex  is  dictated  partly  by  the  conditions 

prevailing  at  its  creation  and  partly  by  the  conditions  surrounding  the  vortex  at  later 

stages  of  its  motion. 

The  results  presented  so  far  have  dealt  with  vortex  pairs  released  at  relatively 

large  depth's  of  submergence  and  have  shown  that  the  computer  code  predicts  results 

which  are  in  good  agreement  with  those  obtained  experimentally  provided  that  a  core 

radius  of  about  r     =  0.07  is  used  in  the  calculations.    It  has  further  been  shown  that 
o 

the  smaller  the  stratification  parameter  the  more  dependent  the  results  become  on  the 
initial  core  readius.  Thus,  the  numerical  analysis  has  not  only  provided  a  power  of 
prediction  for  the  motion  of  the  vortices  but  also  helped  to  delineate  the  proper  value 
of  the  initial  core  radius. 

Having  assertained  the  characteristics  of  the  model  and  the  features  of  the 
velocity,  vorticity,  and  density  variations  it  was  deemed  necessary  to  explore  problems 
of  greater  practicle  and  operational  concern.  It  is  because  of  this  reason  that  the  depth 
at  which  the  vortices  are  generated  was  chosen  to  be  Y/b  =  _  2.0  and  the 
subsequent  motion  of  the  vortices  was  explored.  Figures  B.45  through  B.52  show  the 
densitv,  velocitv,  densitv  fluctuations  and  vorticitv  contours  at  suitable  time  intervals. 
As  the  vortices  approach  the  free  surface  and  begin  to  move  sideways  under  the 
influence  of  their  images,  the  constant  density  contours  begin  to  intersect  with  the  free 
surface  {see  e.g.,  Figure  B.49a)  and  move  sideways  also.  The  intersection  of  these 
density  contours  is  expected  to  give  rise  to  internal-wave-generated  surface  scars.  The 
prediction  of  the  characteristics  of  these  scars  in  terms  of  the  vortices  is  of  importance 
and  requires  that  the  free  surface  condition  be  changed.  This  change  is  to  allow  for 
free  surface  deformation  under  the  influence  of  the  tangential  as  well  as  normal 
velocities  induced  at  the  free  surface  by  the  vortices. 

Additional  calculations  are  underway  to  examine  the  effects  of  free  surface 
deformation,  viscous  diffusion  and  ambient  turbulence  for  various  degrees  of  initial 
stratification.  The  results  presented  herein  and  the  initial  comparisons  are  extremely 
encouraging  and  are  expected  to  lead  to  a  better  understanding  of  this  extremely 
complex  and  challenging  phenomenon. 
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V.  CONCLUSIONS 

The  investigation  reported  herein  warranted  the  following  conclusions: 

(1)  The  motion  of  trailing  vortices  generated  by  submerged  bodies  in  a  stratified 
or  unstratified  medium  can  be  effectively  modeled  by  the  developed  computer 
algorithm: 

(2)  The  computer  model  has  aided  in  the  delineation  of  the  proper  value  for  the 
initial  vortex  core  radius; 

(3)  The  initial  vortex  core  radius  and  vorticity  strength  are  the  primary  factors 
governing  the  stability  of  the  numerical  results; 

(4)  The  migration  of  vortices  is  inhibited  primarily  due  to  the  reduction  of  the 
vorticitv  of  the  initial  vortices  and  the  creation  of  countersigned  vorticitv; 

(5)  The  computer  model  could  be  made  more  realistic  through  the  inclusion  of  a 
deformable  free  surface  in  the  computations. 
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APPENDIX  A 
VORTEX  MOTION  COMPUTER  MODEL 

c 

C  PROGRAM  VORTEX 

C 

C 

C       20  NOVEMBER  86  PROGRAMMER  :  BRIAN  S.  MILLER 

C 

C***************** ************************** ************************** 

c 

C     PURPOSE : 

C 

C     THIS  PROGRAM  NUMERICALLY  SIMULATES  THE  RISE  OF  A  TRAILING 

C      VORTEX  PAIR  SHED  OFF  A  LIFTING  SURFACE  SUBMERGED  IN  A 

C      DENSITY  STRATIFIED  MEDIUM. 

C --- - - 

c 

C  TO  RUN  THIS  PROGRAM: 

C 

C  1)  GO  TO  "INPUT  PARAMETERS"  SECTION  OF  CODE  AND  CHANGE  ANY  PROBLEM 

C  PARAMETERS  (  SP,  FV,  DT ,  ZTMAX,  RO ,  MP,  NP ,  IS AVE)  AS  DESIRED, 

C 

C  .2)  SET  MAXIMUM  TIME  IN  DATA  STATEMENT  MAXSTP; 

C  EX:   DATA  MAXSTP  /N/ 

C  (  WHERE  N=  MAXIMUM  NUMBER  OF  TIME  STEPS) 

C 

C  3)  SET  INTERVALS  FOR  PLOTTING  BY  ENTRIES  IN  DATA  STATEMENT  NPLOT ; 

C  EX:   DATA  NPLOT/  1,2,8,12,20,40,50,75/ 

C  (  THIS  PLOTS  OUTPUT  AT  THE  1ST , 2ND , 8TH , 12TH, 40TH, 50TH , AND 

C  7  5TH  PRQGRAM  TIME  STEP) 

C 

C  4)  CHOOSE   WHAT  OUTPUT  DEVICE  TO  USE,  EITHER  TEK  618  OR  COHPRES 
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C  (METAFILE  OUTPUT)  BY  GOING'  TO   GRAPHIC  DEVICE  SECTION  OF  PROGRAM 

C  AND  COMMENTING  OUT  WHICH  OF  TWO  DEVICES  YOU  DONT  WANT  TO  USE. 

C 

C    5)  COMPILE  PROGRAM  UNDER  VS  FORTRAN  USING  THIS  COMMAND; 

C  FORTVS   VORTEX    (OPT(3) 

C 

C  (  NOTE:  THIS  OPTIMIZES  CODE  AT  LEVEL  3  ) 

C 

C    6)  GO  TO  APPROPRIATE  TERMINAL. AND  RUN  UNDER  DISSPLA. 

C  A.   WHEN  THE  DISSPLA  EXEC  PROMPTS  YOU  ASK  FOR  5  CYLINDERS  OF 

C  TEMPORARY  DISK  STORAGE  AND  ALSO  YOU  MAY  WISH  TO  ISSUE 

C  B.   WHEN  THE  EXEC  PROMPTS  YOU  YOU  MAY  WISH  TO  ISSUE: 

C  FILEDEF   06   DISK  VOR10UT  LISTING  (PERM 

C  IN  ORDER  TO  REROUTE  THE  PRINTED  OUTPUT  TO  A  DISK  FILE 

C  NAMED:  VOR10UT  LISTING  Al  . 

C  C.   PROGRAM  WILL  NOW  RUN  SENDING  OUTPUT  TO  YOUR  DESIGNATED 

C  DEVICES. 

C-- —  - 

C 

C  ADDITIONAL  PROGRAM  NOTES: 

C 

C  1)  THIS  PROGRAM  COMPUTES  DENSITY  AND  VORTICITY  IN  A  (M-2  X  N-2)  GRID 

C  BY  SOLVING  THE  BOUYANTLY  SCALED  FORM  OF  THE  NAVIER-STOKES 

C  EQUATIONS  VIA  AN  "UPWIND  DIFFERENCING"  FINITE  DIFFERENCE  METHOD. 

C   .  •  .  .  '  •       "•'•'.•'■ 

C  2)  THE  VELOCITY  FIELD  IS  COMPUTED  USING  A  FORM  OF  THE  BIOT-SAVART 

C  LAW  AT  EACH  NODAL  POINT  INCLUDING  THE  CONTRIBUTION  OF  THE 

C  "IMAGE  VORTICIES  "  ADDED  TO  INCLUDE  THE  FREE  SURFACE  EFFECT. 

C 

C  3)   THIS  VERSION  OF  THE  PROGRAM  INCLUDES  A  DISSPLA  ROUTINE 

C  TO  PLOT  THE  DENSITY  PERTUBATIONS  , CONSTANT  DENSITY, 

C  VORTICITY,  AND  VELOCITY  FIELD  AT  EACH  SELECTED  TIME  STEP. 

C  THESE  SELECTIONS  FOR  PLOTTING  ARE  INDICATED  BY  ENTRIES 

C  IN  THE  ARRAY  NSTEP.  MAXIMUM  NUMBER  OF  TIME  STEPS  IS  INDICATED 

C  BY  AN  ENTRY  IN  THE  ARRAY  MAXSTEP . 
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c 
c 

C  4)   THE  COMPUTATIONAL  GRID   (M-2)X(N-2)   IS  USED  FOR  ALL  COMPUTATIONS  ' 

C  EXCEPT  FOR  VELOCITY  CALLS  TO  THE  3I0T-SAVART  SUBROUTINES. 

C  TWO  ADDITIONAL  GRIDS,  ONE  FOR  U  AND  ONE  FOR  V,  ARE  SHIFTED  H/2 

C  UNITS  IN  THE  X  AND  Y  DIRECTION  RESPECTIVELY.  THE  SUBROUTINES 

C  COMPUTE  U  AND  V  VELOCITIES  ON  THESE  GRIDS  USING  VORTICITY  FROM 

C  COMPUTATIONAL  GRID.  THE  VELOCITIES  ARE  THEN  AVERAGED  ACROSS  THE 

C  COMPUTATIONAL  GRID  POINTS  TO  BE  USED  FOR  THE  FINITE  DIFFERENCING. 

C 

C  5)   THIS  VERSION  OF  THIS  PROGRAM  WILL  BE  WRITTEN  IN.  SINGLE  PRECISION 

C  AND  IN  AN  OPTIMIZED  FORM  OF  THE  ORIGINAL  CODE. 

C 

C  6)   THIS  PROGRAM  HAS  THE  OPTION  OF  SAVING  THE  RHO  AND  ZETA  ARRAYS 

C  AND  OTHER  PROGRAM  CONSTANTS  IN  ORDER  TO  RESTART  THIS  SIMULATION 

C  AT  THE  TIME  A  PREVIOUS  RUN  HAD  FINISHED.  PROGRAM  VORRS  IS 

C  SPECIFICALLY  WRITTEN  TO  INPUT  THE  DATAFILE  WRITTEN  BY  THIS 

C  PROGRAM  AND  RESTART  THE  SIMULATION  WITH  THE  SAME  PARAMETERS, 

C 

C  TO  ACCESS  THIS  SAVE  FEATURE  SET  VARIABLE  ISAVE  =1  .  A  FILEDEF 

C  OF  FORM  . . .  FILEDEF  08  DISK  <FN>  <FT>  <FM>  . . .  MUST  BE  ISSUED 

C  AT  THE  APPROPRIATE  TIME  IN  EXECUTION. 

C  (NOTE:  TO  DISABLE  THIS  SAVE  FEATURE  SET  ISAVE  =0  ) 

C 

C  •'■    .     ■•..,••:.'■•     :-.:■•.':..•;           .■.'■.. 

C  7)   TO  CHANGE  THE  MESH  DENSITY  OF  THE  COMPUTATIONAL.  DOMAIN  THE 

C  USER  SHOULD  DO  THE  FOLLOWING: 

C  A.   ADJUST  M  AND  M 

C  B.   ADJUST  INITIAL  VORTEX  POSITION:  MP,  MP 

C  C.   ADJUST  THE  DIMENSIONS  OF  THE  W  ARRAY  IN  SUBROUTINES 

C  CPLOT,  ZPLOT,  AND  DRPLOT  TO  BE  (MW,NW) 

C  D.   IF  THE  GRID  IS  LARGER  THAN  50  X  100  THEN  ADJUST  THE 

C  DIMENSION  STATEMENTS  IN  THE  MAIN  PROGRAM  AND  IN 

C  SUBROUTINES  ZPL,  VPLOT  AND  V2PLOT  TO  THE  NEW  GRID 

C  DIMENSION. 
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c 
c 
c 

c 

C      MAJOR  VARIABLES  USED  IN  THIS  PROGRAM; 

C 

C  REAL  VARIABLES; 

C  T     =  NONDIMENSIOMAL  TIME 

C  DT    =  TIME  STEP  • 

C  H     =  NONDIMENSIONAL  GRID  LENGTH  PARAMETER=  DX  =  DY 

C  RO      NONDIMENSIONAL  VORTEX  CORE  RADIUS 

C  XWIDTH=  NONDIMENSIONAL  WIDTH  OF  COMPUTATIONAL  AREA 

C  ZTMAX  =  MAXIMUM  VORTEX  STRENGTH 

C  SP      NONDIMENSIONAL  STRATIFICATION  PARAMETER 

C  FV      NONDIMENSIONAL  FROUDE  NUMBER 

C  BVND   =  NONDIMENSIONAL  BRUNT-VALSA  FREQUENCY 

C 

C  INTEGER  VARIABLES; 

CM     =  NUMBER  OF  NODES  IN  X  DIRECTION 

C  N     =  NUMBER  OF  NODES  IN  Y  DIRECTION 

C  MP    =  X  NODAL  POINT  OF  VORTEX  START  POSITION 

C  MP    =  Y  NODAL  POINT  OF  VORTEX  START  POSITION. 

C  MM     =  X  ARRAY  DIMENSION 

C  NN    =  Y  ARRAY  DIMENSION 

C  NFIG  =  PLOT  NUMBER 

C  NSTEP  =  ITERATION  NUMBER 

C  \ 

C  ARRAY  VARIABLES; 

C  RHO    =  DENSITY  PERTUBATIONS  AT  EACH  NODE  DUE  TO  DYNAMICS 

C  RHOI   =  INITAL  VALUES  OF  DENSITY  AT  EACH  NODE  DUE  TO  GRADIENT 

C  RHOT   =  SUM  OF  BOTH  DENSITY  EFFECTS  AT  EACH  NODE 

C  RHNEW  =  NEW  VALUE  OF  RHO  AT  NEXT  TIME  STEP 

C  ZETA   =  VORTICITY  AT  EACH  NODE 

C  ZTMEW  =  MEW  VALUE  OF  VORTICITY  AT  NEXT  TIME  STEP 
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C  U     =  X  COMPONENT  OF  VELOCITY  AT  EACH  NODAL  POINT 

C  V     =  Y  COMPONENT  OF  VELOCITY  AT  EACH  NODAL  POINT 

C  U2    =  U  VELOCITYS  COMPUTED  AT  BIOT-SAVART  NODAL  POINTS 

C  V2    =  V  VELOCITYS  COMPUTED  AT  BIOT-SAVART  NODAL  POINTS 

C  X     =  X  COORDINATES  AT  EACH  NODAL  POINT 

C  Y     =  Y  COORDINATES  AT  EACH  MODAL  POINT 

C  XU    =  X  COORDINATE  AT  EACH  U  3.  S.  MODAL  POINT 

C  YU    =  Y  COORDINATE  AT  EACH  U  B.  S.  NODAL  POINT 

C  XV    =  X  COORDINATE  AT  EACH  V  B.  S.  NODAL  POINT 

C  YV      Y  COORDINATE  AT  EACH  V  B.  S.  NODAL  POINT 

C 

C - - 

C 

C  SUBROUTINES  USED  IN  THIS  PROGRAM: 

C 

C  SUBROUTINE  COUT   -  OUTPUTS  NODE  COORDINATES 

C 

C  SUBROUTINE  OUTPT  -  OUTPUTS  VARIOUS  PROGRAM  ARRAYS  WHEN  CALLED 

C 

C  SUBROUTINE  CPLOT  -  CREATES  ARRAY  OF  DENSITY  AT  EACH  POINT 

C  AND  CALLS  DISSPLA  CONTOURING  ROUTINE  TO  PLOT. 

C 

C  SUBROUTINE  ZPLOT  -  CREATES  ARRAY  OF   VORTICITY  AT  EACH  POINT  AND 

C  CALLS  DISSPLA  COUTOURING  ROUTINE  TO  PLOT. 

C  .    '  .-■  .'•  •  '"'-.;  •■.  .-■  ....  .''   ■'-;■'■'  '.  •  ..    :■"'■ 

C  SUBROUTINE  DRPLOT-  CREATES  ARRAY  OF  DENSITY  PERTUBATIONS  AT  EACH 

C  POINT  AND  CALLS  DISSPLA  CONTOURING  ROUTINE 

C  TO  PLOT. 

C 

C  SUBROUTINE  ZPL    -  OUTPUTS  VORTICITY  CONTOURS  AS  A  PRINTER  PLOT. 

C 

C  SUBROUTINE  VPLOT  -  CREATES  COMPLEX  VELOCITY  VECTORS  AT  EACH 

C  NODAL  POINT  AMD  THEN  CALLS  THE  FOLLOWING 

C  SUBROUTINES  TO  PLOT  THEM  ON  DISSPLA: 

C  A)      V2PLOT 
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C  B)  AMIN 

C  C)  AMAX 

C 

C  SUBROUTINE    COM      -    CALLED   BY   ZPL   TO   CONTOUR  VORTICITY   ON   PRINTER. 

C 

C 

C      FUNCTIONS    USED    IN   THIS   PROGRAM: 

C 

C  FUNCTION  VBS    -    COMPUTES   THE   V   COMPONENT   OF   VELOCITY   USING 

C  THE   BIOT-SAVART   LAW. 

C 

C  FUNCTION   UBS    -      COMPUTES   THE   U   COMPONENT   OF   VELOCITY   USING 

C  THE   BIOT-    SAVART   LAW. 

C 

C  FUNCTION   PNEW   -   TAKES   THE    FINITE   DIFFERENCE   TO   STEP   THE 

C  EQUATIONS    IN  TIME. 

C 

£******** *x ****** **x *********** *A^ 

DIMENSION   RHOI(50,100) ,U(50,100) ,V(50,10O) ,ZETA(50 , 100) 

*  ,ZTNEW(50,100),RHNEW(50,100),NPLOT(8),X(50),Y(100),RHO(50,100) 

*  ,XU(50),YU(100) ,XV( 50) ,YV( 100) ,U2( 50, 100) ,V2(50,100) 
CHARACTER'S   LABEL 

COMMON  WK( 15000) 

COMMON    /PM/    DT,I,J,UF/UFA,UB,UBA,VF/VFA,VB/VBA,VPT/FV 

COMMON    /PARM/    H,M,N,PI 

DATA  MAXSTP/40    / 

DATA  NPLOT/ 1,10 ,20 ,30, 40, 50, 60, 70/ 

C       -  \ 

C (    CALL   GRAPHIC   OUTPUT   DEVICE    ) 

CALL    COMPRS 
C                CALL   TEK613 
C - - - 

c 

MM=    50 
MN=    100 


c 

C  INPUT  PROGRAM  PARAMETERS  

C 

DT=  2.000E0 
T=0.0 

XWIDTH  =5.0 

R0=  0.09E0 

S?=  1.20 

FV=0. 029765 

ZTMIN=  1.0E-3 

M=  44 

N=  84 

ISAVE=  1 
C DESIGNATE  NODAL  POINTS  TO  PLACE  INITAL  VORTEX  AT  (MP,NP) 

MP=6 

NP=66 

C COMPUTE  ADDITIONAL  PARAMETERS -- 

C 

C...j.  DENSITY  GRADIENT  IN  Y  DIRECTION  INITALLY  =  DRHO 

DRHOI=-SP*SP*FV*FV 
C NONDIMENSIONAL  B.  V.  FREQUENCY  =  BVND 

BVND=  SP*SP*FV*FV 
C MAXIMUM  VORTICITY  AT  ANY  POINT  =  ZTMAX 

ZTMAX=-(FV/(R0*R0)) 
C GEOMETRIC  CONSTANT  PI 

PI=  4.0E0*ATAN(1.0E0)   ■ 
C GRID  LENGTH  PARAMETER  =  H 

H=  XWIDTH/ (M-4) 
C 

C COMPUTE  MISC  M  AMD  N  PARAMETERS 

MM1=  M-l 
MP1=  M+l 
MM 2=  M-2 
NP1=N+1 
NM1=N-1 
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NM2=  N-2 
MW  =  (2*M)  -  7 
NW  =  NM2  -  1 
C 

C INITALIZE  COORDINATE  SYSTEM  ORIGINS  ■ 

C 

X(l)=  -H 
Y(l)=  H 
XU(1)=  -1.5E0*H 
YU(1)=  H 
XV(1)=  -H 
YV(1)=   1.5E0*H 
C 

C COMPUTE  COORDINATES 

C   (  NOTE:  COMPUTATIONAL  GRID  IS  IN  4TH  QUADRANT  +X  -Y  ) 
C 

DO  10  1=  2,MP1 

XU(I)=  XU(I-1)+H 
XV(I)=XV(I-1)+  H 

10  X(I)  =  X(I-l)  +  H 
DO  11  J  =  2,NP1 

YU(J)  =  YU(J-l)  -  H 
YV(J)  =  YV(J-l)  -  H 

11  Y(J)  =  Y(J-l)  -  H 
YHT  =  ABS(Y(N-2)) 

C OUTPUT  INITAL  PARAMETERS  OF  THIS  RUN  OF  PROGRAM. 

C 

WRITE(6,1275) 

WRITE(6,1500) 

WRITE(6,1535)  DT,H 

WRITE (6, 1540)  XWIDTH,YHT 

WRITE(6,1550)  SP,FV 

WRITE (6, 1560)  M,N 

WRITE (6, 1570)  BVND 
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C OUTPUT  COORDINATES  OF  NODES 

WRITE(6,1580)  X(MP),Y(NP) 

WRITE(6,1590) 

CALL  COUT(X,Y,MM,NN) 
C 

c  SPECIFY  INITAL  VORTICITY  AMD  DENSITY  FIELDS  

C 

C    (  FOR  GRID  SYSTEM  X,Y   COMPUTATIONAL  GRID   IN  THE  4TH  QUADRANT  ) 

C     (  TERMS  ARG1  -  ARG4  ARE  IN  QUADRANTS  1  -  '4  ) 


C 


ARGMN=  30.0 
XCENT=  X(MP) 
YCENT=  Y(MP) 
DEN=  2.0E0*R0*R0 
ZETA(MP,NP)=  ZTMAX 
DO  30  J=1,N 

DO  30  1  =  1, M 
TRM1=0.0E0 
TRM2=0.0E0 
TRM3=0.0E0 
TRM4=0.0E0 

ARG1=  ((X(I)-XCENT)**2  +  (Y( J)+YCENT)**2)/DEN 

IF  (ARG1  .GT.  ARGMN)  GO  TO  12 
TRM1=  EXP(-ARGl) 
12     ARG2=  ((X(I)+XCENT)**2  +  (Y( J)+YCENT)**2)/DEN 
IF(ARG2  .GT.  ARGMN)  GO  TO  14 
TRM2=  EXP(-ARG2) 
14     ARG3=  ((X(I)+XCENT)**2  +  (Y( J) -YCENT)**2)/DEN 
IF(ARG3  .GT.  ARGMN)  GO  TO  16 
TRM3=  EXP(-ARGS) 
16     ARG4=  ((X(I)-XCENT)**2  +  ( Y( J) -YCENT) **2 ) /DEN 
IF(ARG4  .GT.  ARGMN)  GO  TO  13 
TRM4=  EXP(-ARG4) 
18      IF((J.EQ.NP)  .AMD.  (I  .EQ.MP))   GO  TO  20 
ZETA(I,J)=ZTMAX*(-TRM1  +TRM2  -TRM3  +TRM4) 


IF(ABS(ZETA(I,J))  .LT.  ZTMIN)   ZETA(I , J)=0 . OEO 
20      RHOI(I,J)=   DRHOI*Y(J) 

RHO(I,J)=  O.OEO 
30      CONTINUE 

WRITE(6,1250) 
C     CALL  OUTPT ( ZETA , ' ZTA ( I , J ) ' , MM , NN , 0 ) 
C     CALL  OUTPT (RHOI, 'R0I(I,J) ' ,MM,NN,0) 
C 

C '  INITALIZE  TIME  STEP  COUNTERS 

NFIG=0 
NSTEP=0 
C 
C**********  ENTRY  POINT  OF  TIME  LOOP  *********** 

60     CONTINUE 
C 

C COMPUTE  U  VELOCITY  FOR  ALL  NODES  IN  U  VELOCITY  GRID 

C 

DO  70  1=  1,MP1 
DO  70  J=  1,N 
70     U2(I,J)=  UBS(XU(I) ,YU(J)7MM,NN,X,Y,ZETA) 
C 

C , COMPUTE  V  VELOCITY  FOR  ALL  NODES  IN  V  VELOCITY  GRID 

C 

DO  100  1=  1;M 

DO  100  J=  1,-NPI 
100    V2(I,J)=  VBS(XV(I) ,YV(J) ,MM,NN,X,Y,ZETA) 
C 

C INTERPOLATE  VELOCITIES  FROM  VELOCITY  GRIDS  TO  NODAL 

C  POINTS  ON  COMPUTATIONAL  GRID 

VMIN=  1.0E-08 
DO  120  1=  1,M 
DO  120  J=  1,N 

U(I,J)=  (U2(I,J)  +  U2(I+1,J))/2.0E0 
V(I,J)=  (V2(I,J)+  V2(I,J+1))/2.0E0 
IF  (ABS(U(I,J))  .LT.  VMIN)   U(I,J)=  O.OEO 
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IF  (ABS(V(I,J))  ,LT.  VMIN)   V(I,J)=  O.OEO 
120     CONTINUE 
C     IF  (NSTE?  .EQ.  0)   CALL  OUTPT(V, ' V( I , J)   ' ,MM,NN,1) 
C     IF  (NSTE?  .EQ.  0)   CALL  OUTPT  (U ,  '  U(  I ,  J)   \MH,NN,1) 
C 
C 

C 

C  INCREASE  TINE  COUNTER  3Y  ONE. 

C    CHECK  ITERATION  COUNTER  AND  TERMINATE  IF  NSTE?  IS  MORE  THAN  MAX 

C    COMPUTE  VORTICITY  AND  DENSITY  AT  ALL  INTERIOR  POINTS 

C         (  STORED  AS  ELEMENTS  OF  ZTNEW  AND  RHNEW  ) 


NSTEP=NSTEP+1 

IF  (NSTEP  .GT.  MAXSTP)  GO  TO  1700 

T=  T+DT 

DO  600  1=2, MM1 

DO  600  J=2,NM1 
VPT=  V(I,J) 

UF=  (U(I+1,J)+U(I,J))/2.0E0 
UFA=ABS(UF) 

UB=  (U(I-1,J)+U(I,J))/2.0E0 
UBA=ABS(UB) 

VF=  (V(I,J-1)+V(I,J))/2.0E0 
VFA=ABS(VF) 

VB=  (V(I,J+1)+V(I,J))/2.0E0 
VBA=ABS(VB) 

ZTNEW ( I , J ) =  PNEW (ZETA , RHO , MM , MM , - 1 . 0E0 , 0 . 0E0 , 0 . 0E0 ) 
600        RHNEW(I,J)=  PNEW(RHO,RHO,MM,NN,0.0E0,0.0E0,BVND) 
C 
C 

C  --- APPLY  BOUNDRY  CONDITIONS  TO  RHO  AND  ZETA  

C 

C ASSIGN  NEW  VALUES  OF  RHO  AND  ZETA  TO  RHO  AND  ZETA  ARRAYS 

C 
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DO   700    I   =   2, MM 
DO   700   J  =   2,NM1 

ZETA(I,J)=   ZTNEW(I,J) 
700  RHO(I,J)=RHNEW(I,J) 

C CALCULATE  RHO  AND  ZETA  VALUES  AT  I=M-2  AND  J=N-2  BOUNDRIES  ... 

C 

DO  300  1=  2,MM2 

DO  300  J=  2,NM2 

ZETA(MM2,J)=  (ZETA(MM1,J)+  ZETA(M-3 , J) ) /2 . 0E0 
ZETA(I,NM2)=  (ZETA(I,NM1)+  ZETA(I ,N-3 ) )/2 . 0E0 
RHO(MM2,J)=  (RH0(MM1,J)+  RHO(M-3 , J) )/2 . 0E0 
RHO(I,NM2)=  (RH0(I,NM1)+  RHO (I ,N-3) )/2 .0E0 

IF(ABS(ZETA(I, J)) . LT .  ZTMIN)   ZETA( I , J)=0 . 0E0 
IF(ABS(RHO(I,J)).LT.  ZTMIN)   RHO(I,J)=  0.0E0 
800     CONTINUE 
C 

C 

C    PLOT  FLOW  PATTERN  WHEN  NSTEP  =  VALUE  SPECIFIED  IN  NPLOT 

C 

C    THEN  RESTART  ALGORITHM  AT  LINE  60  FOR  A  NEW  TIME  STEP 

C .  .' . 

DO  1100  K=  1,8 

IF  (NSTEP-NPLOT(K))   1100,   1200,   1100 
1100  CONTINUE      .'•■'•■;:        -'■■■:.-, 

GO  TO  60 
1200  NFIG=NFIG+1 
TM=  T*FV 

IF  (NSTE?  .EQ.  MAXSTP  )  CALL  OUTPT(ZETA, ' ZTA( I , J) ' , MM , MM , NSTEP) 
IF  (NSTEP  .EQ.  MAXSTP)  CALL  OUTPT  (RHO,  .' RHO (I , J) J  , MM, NN, NSTEP) 
WRITE (6, 1300) NSTEP, T,TM 

CALL  CPLOT ( RHO , RHOI , MM , MM , NSTEP , TM , MW , NW) 
CALL  VPLOT ( U , V , X , Y , MM , NN , TM ) 
CALL  DRPLOT ( RHO , MM , NN , NSTEP , TM , MW , MW ) 
CALL   ZPL ( ZETA , MM , NN , XWIDTH , YHT , ITER ) 
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CALL  ZPLOT(ZETA,MM,NN,TM,MW,NW) 
1250  FORMAT (//  10X,'INITAL  ZETA  AND  RHO  VALUES  :  ') 
1275  FORMATCl1  ) 

1300  FORMAT (/10X, 'ITERATION  NUH: » , 16, 4X, 'TIME=' ,F9 .3,4X, 'T*  =',F10.6) 
1400  WRITE(6,1500) 
1500  FORMAT (  ///) 

1530  FORMAT ('  \5X,'INITAL  PROGRAM  PARAMETERS:1  ) 

1535  FORMAT ( '0' ,5X, 'TIME  STE?=  ' , F 3 . 4 , 5X, ' GRID  SPACING  H  =  ' ,F8.4) 
1540  FORMAT('0' ,5X, ' NON  DIMENSIONAL  CELL  DIMENSIONS  :  XDIRECTICN  =', 

*  F3.4,5X, 'YDIRECTION  =',F8.4) 

1550  FORMAT ( '0' ,5X, 'STRATIFICATION  PARAMETER  ='/F8.4,5X/ 

*  'FROUDE  NUMBER  =',F8.4) 

1560  FORMAT( '0' ,5X, 'NUMBER  OF  NODES  IN  X  DIRECTION  =',I4,5X, 

*  'IN  Y  DIRECTION  =' ,14) 

1570  FORMAT ( '0' ,5X, ' NONDIMENSIONAL  B.  V.  FREQUENCY  =',F15.9) 

1580  FORMAT ( '0' ,5X, 'X  COORDINATE  OF  VORTEX  START  POSITION  =',F8.4,5X, 

*  'Y  COORDINATE  =' ,F8.4) 

1590  FORMAT (' 01 ,5X, 'OUTPUT  COORDINATES  OF  ALL  OTHER  NODES:1) 
C 

C        LOOP  BACK  TO  STARTING  POINT  OF  ITERATION 
GO  TO  60 

1700  CALL  DONEPL 
C1700  CONTINUE 
C 
C SAVE  ARRAYS  AND  TIME  DATA  TO  RESTART  PROGRAM  LATER - 


C 


IF  (ISAVE  .EQ.  0  )  GO  TO  1900 

WRITE (3,*)  SP,FV,DT,T 

WRITE (3,*)  M,N,MP,NP 

WRITE (3,*)  BVND,H,ZTMAX 

DO  1500  1=  1,M 

DO  1300  J=  1,N 
1800     WRITE(8,*)  ZETA(I,J),  RHO(I,J) 
1900   STOP 

END 
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c 

C  +++++++H 

c 

C  SUBROUTINES : 

C 

C 
C===================================================================== 

REAL  FUNCTION   PNEW(P,Q,MM,NN,A,B,C) 
C 

C  THIS  FUNCTION  INTEGRATES  THIS  EQUATION  WITH  RESPECT  TO  TIME 

C  $$$   DP/DT=  (-D(UP)/DX  -D(VP)/DY  +  A*DQ/DX  +  B*DEL(P)  +C*V  )     $$$ 

C     UTILIZING  AN  UPWIND  DIFFERENCING  SCHEME  ' 

C 

DIMENSION  P(MM,NN) ,Q(MM,NN) 

COMMON  /PN/  DT,I, J,UF, UFA, UB ,UBA , VF , VFA, VB , VBA, VPT , FV 

COMMON  /PARM/  H,M,N,PI 

Pl=  (UF-UFA)*P(I+1,J)  +(UF+UFA-UB+UBA)*P(I,J)  - (UB+UBA)*P(I-1 , J) 

P2=  (VF-VFA)*P(I,J-1)  +(VF+VFA-VB+VBA)*P(I,J)  - (VB+VBA)*P(I , J+l ) 

P3=  Q(I+1,J)  -Q(I~lrJ)      . 

P4=  P(I+1,J)  +P(I-1,J)+P(I,J-1)+P(I,J+1)-4.0E0*P(I,J) 

P5=  2.0E0*H*VPT 

TEMP1=(DT/(2.0E0*H))*((-P1-P2)  +A*P3+  (2  .  0E0/H)*B*P4+C*?5  ) 

PMEW=P(I,J)+TEMP1 
■  • RETURN 

END 
C 

c 

REAL  FUNCTION " UBS ( XU , YU , MM , MN , X , Y , 2ETA ) 
C 

C     THIS  SUBROUTINE  FINDS  BOUNDRY  VALUES  OF  VELOCITY   U 
C  UTILIZING  THE  BIOT-SAVART  LAW 

C 
C   (NOTE:  TERMS  IN  EQUATION  ARE  EVALUATED  CONSECUTIVELY  IN  EACH 
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C         QUADRANT  FROM  ONE  TO  FOUR   ) 

C-- 

DIMENSION  X(MM) , Y(NN) , ZETA(MM,NN) 

COMMON  /PARM/  H,M,N,PI 
C 

UMIN=  1.0E-08 

Ul=  0.0E0 

U2=  0.0E0 

U3=  O.OEO 

U4=  O.OEO 
C 

SIGNP=1.0EO 

SIGNM=-1.0E0 

DA=  H*H 
C 
C  COMPUTE  U  COMPONENT  OF  VELOCITY  USING  THE  B-S  LAW 


C 


DO  20  1=  1,M 

DO  20  J=  1,N 
.—  QUADRANT  ONE  -- 

RSQ1=((XU-X(I))**2  +  (YU+Y(J))**2) 

Ul=  U1+  ((-Y(J)-YU)/(2.0E0*PI*RSQ1))*SIGNM*ZETA(I,J)*DA 
--  QUADRANT  TWO  -- 

RSQ2=  ((XU+X(I))**2  +(YU+Y(J))**2) 

U2=  U2+  ((-Y(J)-YU)/(2.0EO*PI*RSQ2))*SIGNP*ZETA(I,J)*DA 
--  QUADRANT  THREE  -- 

RSQ3=  ((XU+X(I))**2  +(YU-Y(J))**2) 

U3=  U3+((Y(J)-YU)/(2.0E0*PI*RSQ3))*SIGNM*ZETA(I,J)*DA 
--  QUADRANT  FOUR  -- 

RSQ4=  ((XU-X(I))**2  +  (YU-Y(J))**2) 
20     U4  =  U4  +  ((Y(J)-YU)/(2.0E0*PI*RSQ4))*SIGNP*ZETA(I,J)*DA 
UTEMP=  U1+U2  +U3+  U4 

IF  (ABS(UTEHP) .LT.UMIN)  UTEMP=0.0E0 
UBS=  UTEMP 
RETURN 
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END 
C 

c 
c=================================================================== 

c 

REAL  FUNCTION  VBS(XV,YVfMM,NN,X,Y,ZETA) 

C 

C.     THIS  SUBROUTINE  FINDS  BOUNDRY  VALUES  OF  VELOCITY  V 

C  UTILIZING  THE  BIOT-SAVART  LAW 

C 

C  (NOTE  :  TERMS  IN  THE  EQUATION  ARE  EVALUATED  CONSECUTIVELY  IN  EACH 

C  QUADRANT  FROM  ONE  TO  FOUR) 


C- 


DIMENSION  X(MM) ,Y(NN) ,ZETA(MM,MN) 

COMMON/PARM/  H,M,N,PI 
C 

VMIN=  1.0E-08 

Vl=  0.0E0 

V2=  0.0E0 

V3=  0.0E0 

V4=  0..0E0 
C 

SIGNP=1.0E0 

SIGNM=-1.0E0 

.  ;-.•:'■  -da=  h*h 

C 

C  COMPUTE  V  COMPONENT  OF  VELOCITY  USING  THE  B-S  LAW 

C  \ 

DO  20  1=  1,M 

DO  20  J=  1,N 
C    --  QUADRANT  ONE  -- 

RSQ1=((XV-X(I))**2  +  (YV+Y(J))**2) 

71=  V1+  ((XV-X(I))/(2.E0*PI*RSQ1))*SIGNM*ZETA(I,J)*DA 
C    --  QUADRANT  TWO  -- 

RSQ2=  ((XV+X(I))**2  +(YV+Y(J))**2) 
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V2=  V2+  ((XV+X(I))/(2.0E0*PI*RSQ2))*SIGNP*ZETA(I,J)*DA 
C.   --  QUADRANT  THREE  -- 

RSQ3=  ((XV+X(I))**2  +(YV-Y(J))**2) 

V3=  V3+ ( ( XV+X ( I ) ) / ( 2 . 0E0*PI*RSQ3 ) ) *SIGNM*ZETA ( I , J ) *DA 
C    --  QUADRANT  FOUR  -- 

RSQ4=  ((XV-X(I))**2  +  (YV-Y(J))**2) 
20     V4  =  V4  +  ((XV-X(I))/(2.0E0*PI*RSQ4))*SIGNP*ZETA(I,J)*DA 
VTEMP  =  V1+V2  +V3+  V4 
IF(ABS(VTEMP) .LT.VMIN)   VTEMP=0.0E0 
VBS=  VTEMP 
RETURN 
END 
C 

c 

c - - 

c 

c 
c 

c 


SUBROUTINE  COUT(X, Y,MM,NN) 


50 


REAL  X(MM),Y(NN),H,PI 

COMMON/PARM/H , M , N , PI 

WRITE (6 ,50) 

WRITE (6,*)  '         ARRAY  X1 

WRITE (6,*) 

WRITE(6,*)  (X(I),  1=1, M) 

WRITE(6,*) 

WRITE (6,*) 

WRITE (6,*)  '         ARRAY  Y 

WRITE  (6,*) 

WRITE(6,*)  (Y(J),  J=1,N) 

WRITE(6,*) 

WRITE(6,*) 

FORMAT (///) 

RETURN 
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END 


SUBROUTINE  OUTPT ( A , LBL , MM , NN , KAL) 
C 

C    THIS  SUBROUTINE  PRINTS  OUT  AN  (NXM)  ARRAY  WITH  A  LABEL 
C  IN  A  RECTANGULAR  (  2 ,M-2)X(2 ,N-2)  GRID 

C 

DIMENSION  A(MM/NN) 

COMMON/PARM/  K , M, N , PI 

CHARACTER*8   LBL 

NM2=  N-2 

MM2=  M-2 

WRITE(6,50) 

WRITE (6,*) 

WRITE (6,*)  '       OUTPUT  OF  ARRAY:    ' , LBL , '    ITER  NUM : ' ,KAL 

WRITE(6,*) 

DO  10  J=  2,NH2 

WRITE (6,*) 

WRITE(6,*)  (A(I,J),  I=2,MM2) 
10   CONTINUE 
50   FORMAT('l'  ) 

RETURN 

END 
C 


SUBROUTINE  CPLOT ( RHO , RHOI , MM , NN , NSTEP , TM , MW , NW ) 

COMMON  WK(15000) 

COMMON  /PARM/  H,M,N,PI 

REAL  W(81,81) 

REAL  RHO(MM,NM) ,RHOI(MM,NN) ,H , PI , RHOT( 50 , 100) 

DATA  ISCALE/4HSCAL/ 
C 

C... CREATE  ARRAY  OF  DENSITY  CHANGES  AT   EACH  NODE  (DRHO) 
C 
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MM2=  M-2 
MM2=  N-2 
DO  10  1=  2,MM2 

DO  10  J=  2,NH2 
JJ=  (NM2+1)-J 
11=  I+(M-5) 

RHOT ( I , J ) =RHO ( I , J ) +RHO I ( I , J ) 
W(II/JJ)=  RHOT(I,J) 
10      CONTINUE 

DO  20  1=  3,MM2 
11=  (M-l)-I 

DO  20  J=  2,NM2 
JJ=  (MM2+1)  -J 
20       W(II,JJ)=  RHOT(I,J) 
C 
C        CALL  OUTPT(RHOT, 'RHT(I, J) ' ,MM,NN,NSTEP) 

C. CALL  SUBROUTINE  CONTR  TO  PLOT  THESE  PERTUBATIONS 

C 

CALL  RESET  (3HALL) 
CALL  PAGE(8.7,11.2) 
CALL  PHYSOR(2.0,3.00) 
CALL  AREA2D  (5.0,6.0) 
CALL  HEIGHT ( .2  ) 
CALL  INTAXS 
■  CALL  XNAME  ('X/BO',4)  . 

CALL  YNAME  ('Y/BO1 ,4) 
CALL  XTICKS(2) 
CALL  YTICKS(2) 

CALL  GRAF  (-5  .0  , 1 .  ,  5  . ,  10  .' ,  1 .  ,0. ) 
CALL  MIXALF( ' INSTRU ' ) 

CALL  MESSAG( ' T(EH . 5 )* (EXHX)  =$',100,1.5,6.5) 
CALL  RE ALNO ( TM , 1 0 6 ,  ' ABUT '  , ' ABUT ' ) 
CALL  FRAME 
CALL  BCOMON  (15000) 
CALL  COMMIM(2.0) 
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CALL  CONANG  (20.) 
CALL  CONMAK  (W,MW,NW, ISCALE) 
CALL  HEIGHT (.08) 

CALL  CONLIN  (0  ,  5HSOLID  ,  '  LABELS  '  ',  2  ,  5  ) 
CALL  CONLIN  ( 1 , 4HDASH, 3HN0LABELS , 1 , 3 ) 
CALL  RASPLN(0.25) 
CALL  COMTUR  (2 , 6HLABELS , 4HDRAW) 
CALL  HEIGHT ( .2) 
CALL  COMPLX 
CALL  RESET ('HEIGHT1 ) 
CALL  RESET ( 'COMPLX1 ) 
CALL  ENDPL  (0) 
RETURN 
END 
C 

c - - ---- 

c 

SUBROUTINE  ZPLOT (ZETA , MM , NN , TM , MW , NW ) 

REAL  ZETA(MM,NN) ,H,PI 

COMMON  WK(15000) 

COMMON  /PARM/  H,M,N,PI 

DIMENSION  W(81,81) 

DATA  ISCALE/4HSCAL/ 
C 
C... CREATE  ARRAY  FOR  PLOTTING  OF  ZETA  IN  ACTUAL  COMPUTATIONAL  AREA, 


C 


MM2=  M-2 

NM2=N-2 

DO  10  1=  2,MM2 
11=  I+(M-5) 
DO  10  J=  2,MM2 
JJ=  (NM2+1)-J 
10       W(II,JJ)=ZETA(I, J) 
DO  20  1=  3,MM2 
11=  (M-l)-I 
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DO  20  J=  2,NM2 

JJ=  (NM2+1)  -  J 
20        W(II,JJ)  =  -ZETA(I,J) 
C 

C CALL  SUBROUTINE  COMTR  TO  PLOT  VORTICITY, 

C 


CALL  RESET  (3HALL) 

CALL  PAGE (8. 7 ,11. 2) 

CALL  PHYSOR(2.0,3.00) 

CALL  AREA2D  (5.0,6.0) 

CALL  HEIGHT (.2  ) 

CALL  INTAXS 

CALL  XNAME  ( ' X/BO ' ,4) 

CALL  YNAME  ('Y/BO1 ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  GRAF  (-5.0,1 . ,5 . ,10. ,1. ,0 . ) 

CALL  MIXALF( 'INSTRU' ) 

CALL  MESSAG('T(EH.5)*(EXHX)  =$',100,1.5,6.5) 

CALL  REALNO(TM,106, 'ABUT' , 'ABUT' ) 

CALL  FRAME 

CALL  BCOMOM  (15000) 

CALL  CONANG  (60. ) 

CALL  COMMIM(2.0) 

CALL  CONMAK  (W,MW,NW, I SCALE) 

CALL  HEIGHT (.03) 

CALL  CONLIN  (0 , 5HSOLID , ' LABELS ' , 2 , 5 ) 

CALL  CONLIN  ( 1 , 4HDASH, ' LABELS ', 1 , 3 ) 

CALL  RASPLN(0.25) 

CALL  CONTUR  (2 , 6HLABELS , 4HDRAW) 

CALL  HEIGHT ( .2) 

CALL  COMPLX 

CALL  RESET ( 'HEIGHT' ) 

CALL  RESET ( 'COMPLX' ) 
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CALL  ENDPL  (0) 

RETURN 

END 


C 
C 


SUBROUTINE  DRPLOT ( RKO , MM , MM , NSTE? , TM , MW , MW) 
COMMON  WK( 15000) 
COMMON  /PARM/  H,M,N;PI 
REAL  W(81;81) 
REAL  RHO(MM,NN) ,H,PI 
DATA  ISCALE/4HSCAL/ 
C 

C... CREATE  ARRAY  OF  DENSITY  CHANGES  AT   EACH  NODE  (DRHO) 
C 

MM2=  M-2 
NM2=  N-2 
DO  10  1=  2,MM2 

DO  10  J=  2,NM2 
JJ=  (NM2+1)-J 
11=  I+(M-5) 
W(II,JJ)=  RHO(I,J) 
10       CONTINUE 

DO  20  1=  3,MM2 
11=  (M-l)-I 

DO  20  J=  2,MM2 
JJ=  (NM2+1)  -J 
20        W(II,JJ)=  RHO(I,J) 
C 

C CALL  SUBROUTINE  CONTR  TO  PLOT  THESE  PERTUBATIONS 

C 

CALL  RESET  (SHALL) 
CALL  PAGE (8. 7, 11. 2) 
CALL  PHYSOR(2.0,3.00) 
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CALL  AREA2D  (5.0,6.0) 

CALL  HEIGHT ( .2  ) 

CALL  IMTAXS 

CALL  XNAME  ( ' X/BO ' ,4) 

CALL  YNAME  ( 'Y/BO' ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  .GRAF  (-5  .0  , 1 .  ,  5  .  ,  10  .  ,  1 .  ,  0  .  ) 

CALL  MIXALF( ' INSTRU ' ) 

CALL  MESSAG('T(EH.5)*(EXHX)  =$',100,1.5,6.5) 

CALL  REALNO(TM,106, 'ABUT' , 'ABUT' ) 

CALL  FRAME 

CALL  BCOMOM  (15000)  • 

CALL  COMMIN(2.0) 

CALL  CONAMG  (20.) 

CALL  CONMAK  (W,MW,NW, I SCALE) 

CALL  HEIGHT ( .08) 

CALL  CONLIN  (0 , 5HS0LID , ' LABELS ' , 2 , 5) 

CALL  CONLIN  (1 , 4HDASH, 8HN0LABELS , 1 , 3 ) 

CALL  RASPLN(0.25) 

CALL  CONTUR  (2 , 6HLABELS , 4HDRAW) 

CALL  HEIGHT (.2) 

CALL  COMPLX 

CALL  RESET ( 'HEIGHT1 ) 

CALL  RESET ( 'COMPLX' ) 

CALL  ENDPL  (0) 

RETURN 

END 


C 

SUBROUTINE  CON (Z2 , XRANGE , YRANGE , MX , H , MM , NN , ZMAX , ZMIN ) 
C 

C      THIS  SUBROUTINE  GENERATES  A  CONTOUR  PLOT  OF  ZETA  IN  A  RECTANGULAR 
C      DOMAIN  OF  SIZE:   YRANGE  *  XRANGE  ON  THE  PRINTER. 
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CHARACTERS  SYMBOL (8)  ,GRAPH(iOO) 
REAL  LX,LY 

DIMENSION  Z2(HM,NN) , VALUE (7) 

DATA  SYMBOL  /  ' X ' , ' A ' ,  ' - ' , '  ' , '  '  ,  ' + ' , ' 1 ' , '  7 
VALUE (1)=  ZMIN 
VALUE (2)=  ZMIN*. 75 
VALUE (3)=  ZMIN*. 5 
VALUE (4)=  0.0 
VALUE (6)=  ZMAX*.75 
VALUE (7)=  ZMAX 
VALUE(5)=ZMAX*.45 
NY=  0 . 8*NX*YRANGE/XRANGE 
DELX=  XRANGE/NX 
DELY=  YRANGE/NY 
WRITE (6 ,50) 

WRITE (6 , *)  ' NX= ' , NX ,  ' NY= ' , NY , ' XRANGE= ' , XRANGE , 
*  '  YRANGE=  '  ,  YRANGE 

WRITE (6,*) 

DO  11  1=  1,NX 
11     GRAPH(I)=  'T' 

WRITE  (-6,  6)  (GRAPH(I),  1=  1,NX) 
DO  5  JSYMBL  =  1,NY 
Y=  (JSYMBL-0.5)*DELY 
J=  1+Y/H 
LY=  Y-(J-1)*H 
HMLY=  H-LY 

DO  4  ISYMBL  =  1,NX 

X=  (ISYM3L-0.5)*DELX 

I=l+X/H 

LX=  X-(I-1)*H 

HMLX=  H-  LX 

Al=  HMLX*HMLY 

A2=  HMLX*LY 

A3=  LX*LY 
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A4=   LX*HMLY 

Z2C=    (A1*Z2(I,J)    +A2*Z2(I,J+1)    +A3*Z2(I+1 , J+l ) 

*  +A4*Z2(I+1,J))/H**2 
C 

C. . .  DETERMINE  IKE  VALUE  OF  NRANGE  BASED  ON  Z2C 

C 

C PRINT  OUT   CHARACTERS  IN  ARRAY  GRAPH  BY  ROW  

C 

DO  2  K=  1,7 

IF  (Z2C-VALUE(K))  1,  1,  2 

1  NRANGE  =  K 

GO  TO  3 

2  CONTINUE 
NRANGE =  8 

3  GRAPH (  ISYMBL)=  SYMBOL (NRANGE) 

4  CONTINUE 

5  WRITE(6,6)  (GRAPH(I),  1=1, NX) 

6  FORMAT( '  ■ ,10X,100A1  ) 
50      FORMAT ('1'  ) 

DO  60  1=  1,NX 
60     GRAPH(I)=  'B' 

WRITE(6,6)  (GRAPH(I)  ,1=1, NX) 

WRITE(6,*) 

WRITE (6,*)  'MAXIMUM  VALUE  OF  NEGATINE  VORTICITY  =  SYMBOL  (X)  = 

*  ,ZMIN 

WRITE (6,*)  '  75%  OF  MAXIMUN  NEGATIVE  VORTICITY  =  SYMBOL  (A)  ' 
WRITE (6,*)  '  50%  OF  MAXIMUM  NEGATIVE  VORTICITY  =  SYMBOL  (-)  ' 
WRITE (6,*) 
WRITE (6,*)  '  MAXIMUM  VALUE  OF  POSITIVE  VORTICITY  =  SYMBOL  (1)  = 

*  ,  ZMAX 

WRITE(6,*)  '  75%  OF  MAXIMUM  POSITIVE  VORTICITY  =  SYMBOL  (+)' 
RETURN 
END 
C - - 

c 
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SUBROUTINE  ZPL (ZETA , MM , NN , XWIDTH , YHT , ITER ) 

COMMON  /PARM/  H,M,N,PI 

DIMENSION  Z2(50,100) ,ZETA(MM,NN) 
C 

C... CREATE  ARRAY  FOR  PLOTTING  OF  ZETA  IN  ACTUAL  COMPUTATIONAL  AREA, 
C 

MM3=  M-3 

MM3=M-3 

ZMIM=  0.0 

ZMAX=  0.0 

DO  10  1=  1,MM3 

DO  10  J=  -1,NM3 
Z2(I,J)=  ZETA(I+1,J+1) 
C.  FIND  THE  MAXIMUM  POSITIVE  VORTICITY  IN  ARRAY  Z2(MM,MN) 


IF(Z2(I,J)  .GT.  ZMAX)   ZMAX=Z2(I,J) 
C.  .  FIND  THE  MAXIMUM  NEGATIVE  VORTICITY  IN  ARRAY  Z2(MM,NN) 
C 

IF(  Z2(I,J)  .LT.  ZMIN)   ZMIN=  Z2(I,J) 
10     CONTINUE 
C 
C 

IF  (ZMAX  .LT.  1.0E-6  )  ZMAX=1.0 

C CALL  SUBROUTINE  CON  TO  PLOT  VORTICITY 

C  -   .•••;,       .;■■  ■ 

CALL  CON(Z2, XWIDTH, YHT ,  100,H, MM ,NN, ZMAX, ZMIN) 

RETURN 

END 
C -- 

SUBROUTINE  VPLOT ( U , V , X , Y , MM , NN , TM) 

COMMON  /PARM/  H,M,N,PI 

COMPLEX  ZV(50,100),  CV(50,100) 

DIMENSION  U(MM,NN),  V(MM,NN),  X(MM) ,  Y(NN) 

MM2=  M-2 

NM2=  M-2 
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MM3  =  M  -  3 

NM3  =  N  -  3 
C 

C FORM  COMPLEX  U,  V,  AND  Z  ARRAYS  FOR  QUADRANTS  3  AND  4 

C 

DO  10  1=  2,MM2 

DO  10  J=  2,NM2 
11=  1-1 

JJ=  (NM2+1)  -J 
ZV(II,JJ)=  CMPLX(X(I),Y(J)) 
10      CV(II,JJ)=  CMPLX(U(I,J),V(I,J)) 

c 

c 

C PLOT  VECTORS  USING  SUBROUTINE  V2PLOT  

CALL  V2PLOT(ZV,CV;MM3,NM3,TM) 
C 

RETURN 
END 


SUBROUTINE  V2PLOT(Z , C,M,N,TM) 
C 

C   *   SUBROUTINE  V2PLOT  TO  PLOT  MXN  COMPLX  VELOCIT  FIELD  C(M,N)  * 
C   *   OVER  A  COMPLX  DOMAIN  Z(M,N)  * 

C   *  * 

COMPLEX  Z(50;100) , C(50,100) , CMAX 
REAL  X(100) ,Y(100) ,YX(100) ,YN(100) ,XX(100) , 
*  XN(100),X2(100) ,Y2(100) 
INTEGER  PMAX,PMIN 
EXTERNAL  AMAX 
EXTERNAL  AMIN 
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DO  20  J  =  1,N 

DO  10  I  =1;M 

X(I)  =ABS(REAL(C(I,J))) 
Y(I)  =ABS(AIHAG(C(I,J))) 
10        CONTINUE 

CALL  AMAX(Y,M,YMAX,PMAX) 
CALL  AMAX(X,M,XMAX,PMAX) 
YX(J)  =  YMAX 
XX(J)  =  XMAX 
CALL  AMIN(Y,M,YMIN,PMIN) 
CALL  AMIN(X,M,XMIN,PMIN) 
YN(J)  =YMIN 
XN(J)  =XMIN 
20 _     CONTINUE 

CALL  AMAX(YX,N,YMAX,PMAX) 
C        CALL  AMIN(YN,N,YMIN,PMIN) 

CALL  AMAX(XX,N,XMAX,PMAX> 
C        CALL  ANIN(XN/N,XNIN,PMIM) 


C 


C 


CMAX  =  CMPLX(XMAX,YNAX) 
UMAX  =SQRT( ( XMAX) **2+( YMAX) **2) 
DO  24  I  =1,M 
DO  22  J=1,N 

C(I,J)  =  C(I,J)/UMAX 
22       CONTINUE 
24     CONTINUE 

DO  40  J  =1,N 

DO  30  I  =1,M 

X(I)  =REAL(Z(I,J)) 
Y(I)  =AIMAG(Z(I,J)) 
30        CONTINUE 

CALL  AMAX(Y,M,YMAX,PMAX) 
CALL  AMAX(X/M,XMAX,PMAX) 
YX(J)  =  YMAX 
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XX(J)  =  XMAX 

CALL  AMIN(Y,M,YMIN,PMIN) 

CALL  AMIN (X,M,XMIN, PMIN) 

YN(J)  =YMIN 

XN(J)  =XMIN 
40      CONTINUE 

CALL  AMAX(YX,N,YMAX,PMAX) 
CALL  AMIN ( YN , N , YMIM , PMIN ) 
CALL  AHAX(XX,N,XMAX,PMAX) 
CALL  AMIN (XN,N,XMIN, PMIN) 

CALL  RESET(3HALL) 

CALL  PAGE (8. 7, 11. 2) 

CALL  PHYSOR(2.0,3.0) 

CALL  AREA2D ( 5 . , 6 . ) 

CALL  HEIGHT (.2) 

CALL  INTAXS 

CALL  XNAME('X/B0' ,4) 

CALL  YNAME( 'Y/BO' ,4) 

CALL  XTICKS(2) 

CALL  YTICKS(2) 

CALL  GRAF (0.0, 0.5, 5. 0,-10. 0,1. 0,0.0) 

CALL  MIXALFC'INSTRU1 ) 

CALL  HESSAG('T(EH.5)*(EXHX)  =$',100,1.5,6.5) 

CALL  RE ALNO  ( TM ,  1 0  6  ,.  '  ABUT  '  ,  '  ABUT  '  ) 

DO  60  J=1,N 

DO  50  1=1, M 

X(I)  =  REAL(Z(I,J)) 
X2(I)  =  X(I)+REAL(C(I,J)) 
Y(I)=  AIHAG(Z(I,J)) 
Y2(I)  =  Y(I)+AIMAG(C(I,J)) 
50         CONTINUE 

CALL  MARKER  (3) 
CALL  5CLPIC( .001) 
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CALL  CURVE (X,Y,M,-1) 
DO  54  I  =1,M 
.  X11=X(I)  ' 
X12=X2(I) 
Y11=Y(I) 
Y12=Y2(I) 

IF(X12.EO.O. .AND.Y12.EO.0)  GO  TO  54 
CALL  RLVEC  (Xll , Yll , X12 , Y12 , 1121 ) 
54         CONTINUE 
60      CONTINUE 

.  CALL  RESET ('HEIGHT1 ) 
CALL  ENDPL(O) 
RETURN 
END 
C 
C 

SUBROUTINE  AMAX(Y,N, YMAX,PMAX)  . 
C 

Q  ******************************************** 

C        *   SUBROUTINE  TO  COMPUTE  THE  LARGE  ELELENENT 
C        *   IN  A  GIVEN  ROW  OR  COLUMN  IN  ARRAY  ' A ' . 

Q  ******************************************** 

c 

C        ***  VARIABLE  DECLARATION  *** 
C        REAL  AMAX 

INTEGER  PMAX 
DIMENSION  Y(N) 
C 

YMAX=Y ( 1 ) 

DO  5100  1=2, N 

IF(Y(I) .LE.YMAX)  GO  TO  5000 
YMAX=Y ( I ) 
PMAX=I 
5000       ,  CONTINUE  " 
5100     CONTINUE 
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C  MAXROW=POS 

RETURN 

END 
C 
C 
C 

SUBROUTINE   AMIN ( Y , M , YMIN , PMIN) 
C 

C  *      SUBROUTINE   TO   COMPUTE   THE    SMALL   ELEMENT 

C  *      IN   A  GIVEN   ROW   OR   COLUMN   IN  ARRAY    'A' . 

Q  it******************************************* 

C 

C  ***   VARIABLE   DECLARATION   *** 

C  REAL  AMIN 

INTEGER  PMIN 

DIMENSION  Y(N) 


YMIN=Y(1) 

DO   6100    1=2, N 

IF(Y(I).GE.YMIN)    GO   TO    6000 

YMIN=Y(I) 

PMIN=I 

6000 

CONTINUE 

6100 

CONTINUE 

RETURN 

END 
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Figure  B.la    Velocity  Field:  SP  =  0.00,  T*  =  0.05952. 
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Figure  B.lb     Vorticity  Field:  SP  =  0.00. 
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T*   -    0'. 59520 
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Figure  B.2a     Velocity  Field:  SP  =  0.00. 
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Figure  B.2b     Vorticity  Field:  SP  =  0.00. 
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Figure  B.3a     Velocity  Field:  SP  =  0.00. 
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Figure  B.3b     Vorticity  Field:  SP  =  0.00. 
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T*  =    1.78560 
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Figure  B.5a     Velocity  Field:  SP  =  0.00. 
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Figure  B.5b    Vorticity  Field:  SP  =  0.00. 
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Figure  B.6a     Velocity  Field:  SP  =  0.00. 
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Figure  B.6b     Vorticity  Field:  SP  =  0.00. 
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Figure  B.7a     Velocity  Field:  SP  =  0.00. 
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Figure  B.7b     Vorticity  Field:  SP  =  0:00. 
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Figure  B.8a     Velocity  Field:  SP  =  0.00. 
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Figure  B.8b     Vorticity  Field:  SP  =  0.00. 
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Figure  B.9a     Velocity  Field:  SP  =  0.00. 
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Figure  B.9b     Vorticity  Field:  SP  -  0.00. 
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Figure  B.lOa     Velocity  Field:  SP  =  0.00. 
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Figure  B.lOb     Vorticity  Field:  SP  =  0.00. 
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Figure  B.lla     Velocity  Field:  SP  =  0.00. 
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Figure  B.12a     Velocity  Field:  SP  =  0.00. 
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Figure  B.12b     Vorticity  Field:  SP  =  0.00. 
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Figure  B.13a     Velocity  Field:  SP  =  0.00. 
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Figure  B.13b     Vorticity  Fieid:  SP  =  0.00. 
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Figure  B.14     Experimental  and  Numerical  Results  for  Various  Initial  Core  Radii. 
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Figure  B.  15b     Velocity  Field:  SP  -    1.00. 
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Figure  B.  15c     Density  Perturbation  Contours:  SP  =   1.00. 
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Figure  B.l5d     Vorticity  Field:  SP  =    1.00. 
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Figure  B.16a     Constant  Density  Contours:  SP  =   1.00. 


90 


I*  -    0.34380 


2.5    3       3.5     4       4.5    5 
/BO 


Figure  B.16b    Velocity  Field:  SP  =   1.00. 
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Figure  B.16c     Density  Perturbation  Contours:  SP  =   1.00. 
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Figure  B.16d     Vorticity  Field:  SP  =    1.00. 
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Figure  B.20c     Density  Perturbation  Contours:  SP  =   1.00. 
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Figure  B.20d     Vorticity  Field:  SP  =   1.00. 
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Figure  B.21c     Density  Perturbation  Contours:  SP  =   1.00. 
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Figure  B.21d    Vorticity  Field:  SP  =   1.00. 
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Figure  B.22a     Constant  Density  Contours:  SP  =   1.00. 
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115 


r  -    2.40660 


O  -i 


O 
CO 
\ 


CM  - 


ro  - 


TJH    - 


in  - 


co  - 


[N.  - 


co- 


rn - 


o 


-5      -4 


■3     -2-10        1 

X/BO 


2        3 


Figure  B.22c     Density  Perturbation  Contours:  SP  =   1.00. 


116 


V  -    2.40660 


o 

CO 


-1       0        1 

X/BO 


Figure  B.22d     Vorticity  Field:  SP  =   1.00. 
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Figure  B.23a     Constant  Density  Contours:  SP  =   1.00. 
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Figure  B.34a     Constant  Density  Contours:  S?  =  0.50. 
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Figure  B.35a     Constant  Density  Contours:  SP  =  0.50.. 
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